A Monte Carlo Approach for Studying Microphases Applied to the Axial 
Next-Nearest-Neighbor Ising and the Ising- Coulomb Models 
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The equilibrium phase behavior of microphase-forming systems is notoriously difficult to obtain 
because of the extended metastability of their modulated phases. In this paper we present a system- 
atic simulation methodology for studying layered microphases and apply the approach to two proto- 
typical lattice-based systems: the three-dimensional axial next-nearest-neighbor Ising (ANNNI) and 
Ising-Coulomb (IC) models. The method involves thermodynamically integrating along a reversible 
path established between a reference system of free spins under an ordering field and the system of 
interest. The resulting free energy calculations unambiguously locate the phase boundaries. The 
simple phases are not observed to play a particularly significant role in the devil's fiowers. With the 
help of generalized order parameters, the paramagnetic-modulated critical transition of the ANNNI 
model is also studied. We confirm the XY universality of the paramagnetic-modulated transition 
and its isotropic nature. Interfacial roughening is found to play at most a small role in the ANNNI 
layered regime. 

PACS numbers: 64.60.Cn, 64.60.F-,05.10.Ln,75.10.-b 



1. INTRODUCTION 

Lattice models are central to statistical mechanics. 
They strip away the complexity due to packing and help 
reveal the influence of non-geometrical factors on both 
equilibrium and non-equilibrium self assembly. The Ising 
model, for instance, offers a singular window on crit- 
ical phenomena and on gas-liquid coexistence^; Flory- 
Huggins's theory of solvated polymers is core to the 
physics of polymers^; and spin glasses are key sources 
of inspiration for the difficult problem of structural glass 
formation^. If a lattice model of a system exists, it is 
often a good strategy to solve it before embarking on a 
study of more elaborate variants. 

Microphase formation is one such phenomenon that 
could benefit from further consideration of lattice-based 
models. The frustration of short-range attraction - or 
sometimes repulsion^ - by a long-range repulsion, irre- 
spective of the physical and chemical nature of these 
interactions, leads to universal spatially modulated pat- 
terns^. Periodic lamellae, cylinders, clusters, etc. are thus 
similarly found in block copolymers'^*, oil- water surfac- 
tant mixtures^iiS,, charged colloidal suspensions-^, and 
numerous magnetic materialsi^iii^. Microphase forma- 
tion has also been hypothesized to play a role in biological 
membrane organization^** and in the formation of stripes 
in certain superconductors'^"^*', though the microscopic 
interpretation is still debated. The spontaneous nature 
of microphase organization allows for these mesoscale 
periodic textures to find technological success as ther- 
moplastic elastomers' and nanostructure templates^'. 
Obtaining a detailed control over microphase morphol- 
ogy remains, however, notoriously difficulliS^. Anneal- 
ing2^, external fields^'*, strain compression^, addition of 
fullerenes^i^, or complex chemical environments^! are of- 
ten necessary to order diblock copolymers, for instance. 



Understanding how to tune and stabilize microphases 
is essential to broadening their material relevance, yet ex- 
perimental systems provide limited microscopic insights. 
A number of continuous space^^Ti^i and lattice^"— mod- 
els have thus been devised for theoretical and simulation 
studies, and some of which have even become textbook 
materia l^^'^^ . Grasping the equilibrium properties of 
these models is necessary to resolve problems surround- 
ing the non-equilibrium assembly of microphases^"—. 
But although the modulated regime is a key feature of 
these models, it has not been accurately characterized in 
any of them. Even for the most schematic formulations, 
the existing theoretical treatments have only offered lim- 
ited assistance. 

Direct computer simulations have also been unable to 
provide reliable equilibrium informational^. Traditional 
simulation methodologies that facilitate ergodic sampling 
of phase space by passing over free energy barriers, no- 
tably parallel tempering and cluster moves, are of lim- 
ited help in microphase-forming systems. Because of 
the dependence of the equilibrium periodicity on tem- 
perature, sampling higher temperatures leaves the sys- 
tem in a modulated phase with the wrong periodicity; 
and because of the lack of simple structural rearrange- 
ments for sampling different modulations, the efficiency 
of cluster moves is limited. We recently introduced a 
free-energy integration method for simulating modulated 
phases that overcomes this hurdle^**. Here, we detail this 
method and apply it to the study of two canonical three- 
dimensional (3D) spin-based systems: the axial next- 
nearest-neighbor Ising (ANNNI) and the Ising-Coulomb 
(IC) models. Both of these models are known to form 
lamellar phases of different periodicities at low temper- 
ature, but their phase structure is still not completely 
understood. The phase information we obtain by simu- 
lation further allows testing of various theoretical predic- 
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FIG. 1: (Top) Snapshot of the ANNNI antiphase (2) at T = 
2.4, K = 0.7 for a 20 X 20 X 40 lattice. Differently shaded beads 
indicate spins up or down. (Bottom) Notation examples of the 
lamellar phases (23) and (2). 

tions. The plan of this paper is to introduce the models 
(Sect. [2]), the simulation methodology (Sect. |3]), and the 
generalized order and critical parameters (Sect.|4]). After 
discussing the results (Sect.[S|), a short conclusion follows. 

2. MODELS 

Before introducing the models, a clarification of the 
nomenclature for describing layered microphases is in or- 
der. Two conventions for characterizing the periodic- 
ity of lamellar phases coexist in the scientific literature. 
The first compactly identifies a phase with a simple wave 
number q — 1/ A (in units of 27r) , where A is the period 
length. The second, a short-hand form {m^n^) intro- 
duced in Ref. Isil . is less compact but provides a more in- 
tuitive description of the layered phase. In this notation, 
integers are used to describe a lamellar phase formed by 
periodic repetition of patterns of j lamellae of width m 
followed by k lamellae of width n (Fig.[T]). For example, 
phase (oo) is the ferromagnetic phase, phase (2) consists 
of two layers of spins up followed by two layers of spins 
down, and phase (23) has a period of 5. Because ther- 
mal fluctuations blur the layer boundaries, the thickness 
of each lamella is generally not an integer but takes an 
average value 

A mj + nk , . 

2" j+k • ^> 

This notation, which can only represent phases of ra- 
tional periodicity, is well suited for the commensurate 
phases that are here observed. 

2.1. ANNNI Model 

The ANNNI model was first introduced to rationalize 
helical magnetic order in certain heavy rare-earth met- 



g^^g36 -38,52 ^ The simple model's description of the exper- 
imentally observed order is only qualitative^'^, but be- 
cause of its surprisingly complex phase behavior, it is 
now canonical for the study of systems with competing 
interactions^!^. Its Hamiltonian on a simple cubic lat- 
tice 

is expressed for spin variables = ±1 coupled through 
a positive constant J. With the Boltzmann constant 
ks, J/ks sets the temperature T scale. Alignment is 
favored for nearest-neighbor pairs but frustrated 

with relative strength k > for z-axial next-nearest- 
neighbor pairs The exact solution of the one- 
dimensional version of the model provides T = phase 
information for all other dimensions^-: ferromagnetic or- 
der is the ground state for k < 1/2, while the layered 
antiphase (2) minimizes the energy for k > 1/2. A 
mean-field description qualitatively captures the higher- 
dimensional, finite-T features of the model^i^: the sys- 
tem is paramagnetic at high T; it is ferromagnetic at low 
T and k; and modulated layered phases form for suffi- 
ciently high These three regimes join together at 
a multicritical Lifshitz point (kl,Tl) whose special crit- 
ical properties have been predicted by theory^"— and 
verified in simulations^2iS£. High-temperature series ex- 
pansions have also been used to study the paramagnetic 
phase and predict its limit of stability^i^. These pre- 
dictions were confirmed by finite-size critical rescaling for 
the paramagnetic-ferromagnetic (PF) transitional^"— 
and by heat capacit y^^'^^ and generalized susceptibilit}*^ 
measurements for the paramagnetic-modulated (PM) 
transition. For k < kl, the PF transition has Ising 
universality^^i^; while for k > kl, the PM transition 
has been argued to have XY universality^S^^sSE^ but di- 
rect simulation verifications are incomplete^ and the 
results of the high-temperature series expansion analy- 
sis are inconclusive^^'''''. The ferromagnetic- modulated 
(FM) transition is predicted by a Landau-Ginzburg treat- 
ment to be first order with q changing discontinuously 
from 0, and to be tangent to the PF and PM transition 
lines at the Lifshitz poini^i. 

A sequence of commensurate (2^3) phases spring from 
the multiphase point at T = and k = 1/2. The 
structure of the branching processes at low T has been 
carefully studied^, and forms the basis for the low- 
temperature series expansion^. For the rest of the 
modulated regime, approximate theoretical treatments, 
such as an approximate mean-field theory with a soliton 
correctional, an effective-field theory^, and the tensor 
product variational approach (TPVA)2^ have been used. 
Monte Carlo simulations have also been carried out in 
this regime^iH, but the hysteresis resulting from the 
high free-energy barriers that separate modulated phases 
from each other limits accurate determinations of the 
phase boundaries from annealing-based approaches'^ 
Avoiding annealing is thus preferable for accurately lo- 



3 



eating transitions within the modulated regime^. It 
is thought that incommensurate phases could lower the 
transition free energy barriers between different commen- 
surate modulated phases on sufficiently large lattices^^, 
but these phases have not been observed thus far. 



2.2. Ising-Coulomb Model 

The Ising-Coulomb (IC) model, in which the near- 
est neighbor ferromagnetic coupling spin is frustrated by 
long-range Coulomb interaction of relative strength Q, 
was first suggested as a model for the stripe phase behav- 
ior of high-temperature superconductors in two dimen- 
sions^SiZl. It was also adopted as a generic coarse-grained 
description of microphase formation in systems with com- 
peting pair interactions in three dimension a^'^i^^'^^ , and 
used to study the effect of dispersion forces on phase 
transitions in ionic systems^S. Although it is based on 
an Ising model, its Hamiltonian 

Hjc = -JJ2s,s,+QJY;^-^^ (3) 

does not allow ferromagnetic ordering for any Q > 0, i.e., 
an infinitesimally small Coulomb frustration is sufficient 
to induce layering^^. But by analogy with a Landau- 
Ginzburg model with frustration^^iZli^iS^, it is expected 
that any screening of the Coulomb interaction would 
move the onset of modulation to a finite Q^i^. Interest- 
ingly, the Q — >■ oo limit recovers the simple-cubic lattice 
restricted primitive model (LRPM) of Dickman and Stell 
at full occupancy^>2^. 

The one-dimensional T — phase sequence is known 
to be made of equal length blocks of alternating ori- 
entation^i. In higher dimensions, though no rigorous 
demonstration exists, layered phases of integer period- 
icity are also expected to be the ground state at low 
Cfi^. In that regime, an approximate mapping to a one- 
dimensional system seems reasonable. For sufficiently 
large Q, two- and three-dimensional periodic structures, 
i.e., "cylinders" and "clusters", minimize the energy; and 
for Q > w 15.33 antiferromagnetic Neel order is ex- 
pected^^. Mean-field treatment o^^i^° and Monte Carlo 
simulations (for Q < 1)^ describe the paramagnetic- 
modulated (PM) transition. Although the mean-field 
results overestimate the transition temperatur o^^'^" , the 
predictions are nonetheless quite similar to the phase 
behavior obtained from simulations'*'^. Because of the 
long-range isotropic Coulomb interaction, the transition 
is "fiuctuation-induced" first order for any < Q < 
^^^ ,88,89 ^ g^j^^j jp.^ Q modulated phases melt at 

Tc(Q) - Tc(0) - Ql/^ where Tc{0) « 4.51 is the 3D Ising 
simple cubic critical point—. For Q > Qn the continu- 
ous paramagnetic-Neel (PN) transition has Ising univer- 
sality, and at high Q, the critical temperature Tc{Q) ~ 
Tc(oo)Q — Tc(0)22, where the trivial linear dependence re- 
sults from the choice of units and Tc(oo) « 0.515 is known 



from LRPM simulation a^^'^° . A triple point connects the 
paramagnetic, modulated cluster, and antiferromagnetic 
Neel phases at {Qn,Tc{Qn)), where only the mean-field 
estimates Qjv — SG/tt « 11.5 and Tc{Qn) = 1-61 are 
known^^. Within the modulated layered regime proper, 
phases spring out at the boundary between neighboring 
low-temperature ground states of integer periodicity^. 
The process is akin to the springing of phases between 
the antiphase and the ferromagnetic phase in the ANNNI 
model. The simulations in the layered regime capture 
the presence of these phases, but the use of a simulated- 
annealing approach in a strongly hysteretic regime is 
likely to bias the estimates for the transition tempera- 
tures'*'' . 



3. METHOD 

Monte Carlo simulations are used for determining the 
absolute free energy of the different modulated phases. 
The thermodynamic integration, the reference systems, 
and the Monte Carlo sampling details are presented in 
this section. 



3.1. Thermodynamic Integration 

The free energy is obtained from Kirkwood thermody- 
namic integration^*'^^, which involves simulating a sys- 
tem with a Hamiltonian that couples a reference system 
Hamiltonian Hq with that of the system of interest Hi 

Hx^il-X)Ho + XHi. (4) 

For a given A, the Helmholtz free energy F\ obeys 

dFx^ T^dZx^/dHxX 
dX Zx dX \ dX ^ ' 

where Zx is the canonical partition function and (• • 
denotes a canonical average under Hx- The difference 
between the free energy of system of interest Fi and that 
of a known reference system Fq at phase point (Tq, kq) is 
thus 

Jo^ \ dX U 

= / {Hi-Ho)xdX. 
Jo 

In order to obtain reliable numerical results, the inte- 
gration path from A = to 1 must be reversible. No 
first order phase transition may take place along it. Our 
choice of reference system, which is key to the approach, 
is detailed in the next subsection. The numerical integra- 
tion is done by simulating the system at discrete A points 
chosen following a Gauss-Lobatto scheme^''. Because of 
a rapid change in the integration curve as A — >■ 1, the lat- 
ter part of the integral uses logarithmically spaced points 
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FIG. 2: (top left) Thermodynamic integration of the ANNNI model at k — 0.7 and T — 2.5 for phase (2) using sinusoidal and 
square fields as reference, (bottom left) Change in the structure factor peak height along the integration path and the free 
energy results for the two different references, (top right) Thermodynamic integration of the IC model at Q = 0.8 and T = 1.06 
for phase (21*). The integration curve from a fluctuating to a constant magnetization system is shown in the inset, (bottom 
right) The structure factor at different wavevectors demonstrates the preservation of the modulation along the integration path. 



that are densely distributed near A = 1 (Fig. [2]). This ad- 
justment is necessary for accurately capturing the z-axis 
translational degree of freedom of the lattice, whose con- 
tribution is particularly important in small systems. 

In principle, one could investigate the T-frustration 
plane point by point, but that would be computation- 
ally wasteful. The data collection is significantly acceler- 
ated by thermally integrating to nearby temperatures Ti 
or frustrations ki (or, equivalently, Qi) using a known 
state point (To, kq) as reference by 

'"^^ dil/T) (7) 



Fi(Ti,kq) Fi(To,ko) 



Ti 



To 



Fi(To,«;i)-Ti(To,Ko)= / {^) dK, (8) 



where 



and 



dFi 
Ok 



dHi 
Ok 



To 



J SiSj 



(9) 



(10) 



In practice, the free energy results are fitted with a poly- 
nomial of degree three or four. The free energy at any 
point within a relatively short interval is then interpo- 
lated from the parameterized function. 



3.2. Reference System 

In order to guarantee a reversible integration path, the 
reference system should reflect the symmetry of the phase 
under investigation. A good reference system should also 
have a Hamiltonian Hq whose partition function Zq and 
free energy To can be obtained analytically or at least 
with high numerical accuracy. For the lamellar phases 
observed on lattices with N = L^LyL^ sites, we propose 
a reference that has decoupled spins under a z-axial pe- 
riodically oscillating field B{z) with amplitude Bq 



N 



Ho = -BoJ2^^B{z,), 



(11) 



4=1 



similarly to the periodic potential wells confining free 
particles used in Ref. 94. It trivially follows that in a 
system with fluctuating magnetization 



To 
NT 



1 



2 cosh 



BoBjz] 
T 



(12) 



The amplitude Bq should be sufficiently strong to pre- 
vent layer melting and changes of layer periodicity as 
the field is turned off, yet sufficiently weak to allow 
sampling of the integrand^. Fortunately, the relatively 
high free-energy barriers between neighboring modulated 
phases make phase transitions along the integration path 
highly unlikely, even if sections of the path are formally 
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metastable. Due to the broken symmetry between the 
different coordinate axis, we can also, without loss of 
generality, similarly lock the lamellae in a specific orien- 
tation when initializing configurations for the IC model. 

The applied field B{z) needs not be the exact equilib- 
rium profile of the modulated layers as long as the inte- 
gration from B(z) can be done reversibly. For instance, 
either square or sinusoidal fields can be used as reference 
states for the study of modulated phases with integer 
periodicity. The free energy results of both approaches 
agree with high accuracy (Fig. [2]). The equivalence also 
holds in the low-temperature regime, where the ground 
state profile is more akin to a square well than to a pure 
sine function^. Because sinusoidal fields are "soft" in the 
interlayer region, which helps averaging the layer fluctu- 
ations, and because they provide a compact and efficient 
way to describe non-integer periodic lamellae, we use 



3.3. Monte Carlo Sampling 



B{z) — sin{2TTqz + (jio), 



(13) 



where a small phase angle 4>q is added to prevent the 
lattice sites from directly overlapping with the zeros of 
the field. 

The IC model, which must remain charge neutral, re- 
quires that the reference partition function Zq be com- 
puted subject to a fixed magnetization constraint. In 
the infinite system limit this correction is negligible, but 
on a finite lattice it may affect transition temperatures. 
For the paramagnetic phase, the reference system Hamil- 
tonian Hq with Bq = results in Fq/N = —Tin 2 for 
an unconstrained system (Eg. [T^ . but the properly con- 
strained reference system instead has 



NT 



N \N/2 



(14) 



where (^2) binomial coefficient. For N — 12^ x 

24 = 3456 spins, the difference between the two results is 
~ O.OOir, which may be significant because of the small 
entropy differences between layered phases. Calculating 
Fq is not, however, as straightforward for Bq ^ 0. One 
has to define a thermodynamic integration path between 
fluctuating and constant magnetization systems 



(15) 



with A going from to 00. In practice, (^3^) = 
rapidly decays to zero with growing A, and therefore in- 
tegrating to a flnite A of order unity is sufficient. The 
zero-magnetization free energy F'^ is then obtained by 
adding the correction from thermodynamic integration 
to Fq from Eq. [T^] (Fig. [2]). We note, however, that even 
in the small IC systems studied here, the free energy cor- 
rections for different modulated phases are very similar 
for a given temperature. The phase transitions are thus 
only imperceptibly affected by the shift. 



We perform constant T Monte Carlo (MC) simulations 
on a cubic lattice under periodic boundary conditions, us- 
ing N = L^LyL^ = 40^ X 240 spins for the ANNNI model 
and N = 12^ x 24 spins for the IC model, unless otherwise 
noted. Ewald summation is used to compute the long- 
range Coulomb interactions in the IC model^^'^^. The 
phases studied have wave numbers q = n/Lz with inte- 
ger n's, which keeps modulations commensurate with the 
lattice. We initialize the modulated phases with a sinu- 
soidally varying spatial probability of the desired period- 
icity. The system relaxes to the equilibrium spin profile 
for a given xy plane 

Sa;y(2:) = V(s«> (16) 

irrespectively of the initialization scheme as long as it has 
the correct periodicity. 

Basic MC sampling consists of single-spin flips for the 
ANNNI model. Spin exchanges, which enforce charge 
neutrality, are used for the IC model. Phase-space ex- 
ploration gains in efficiency by complementing the basic 
sampling with iterations that take advantage of phase 
symmetry. 

• For the modulated phases, layer swaps allow for 
the individual layer thickness to fluctuate while pre- 
serving the overall periodicity. Multiple layer swaps 
are necessary to alter the periodicity and therefore 
even neighboring modulated phases are well sepa- 
rated in configuration space. 

• Near the PM transition, an anisotropic cluster algo- 
rithm for the ANNNI model^^ and a modified Wolff 
algorithm that considers the corrections from long- 
range interaction for the IC model— are used, in 
order to capture the strong fluctuations. 

• For systems with an applied external magnetic 
field, lattice drifts with respect to the field along 
the z axis help sample the translational degrees of 
freedom. 

For reference point integrations, up to 10^ MC moves 
{N attempted spin flips or exchanges per move) are per- 
formed after 5 x 10^ MC moves of preliminary equilibra- 
tion. For the thermal and frustration integrations, only 
10'* MC moves are necessary, because the free energy is 
not as sensitive to the accuracy of the integration slope as 
it is to its starting point, over the small T intervals con- 
sidered. In the vicinity of the critical transitions, we also 
use the multiple histogram algorithm, in order to obtain 
high precision results with a minimal amount of compu- 
tations^^. The method relies on reweighing the sampled 
configurations at a fixed temperature Tq , typically nearby 
Tc, by the Boltzmann factor difference e-(i/r-i/To)B 
results at neighboring temperatures T— . Our implemen- 
tation uses a logarithmic summation scale, in order to 
avoid sum overflow in large systems^^. 
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FIG. 3: Optimized and un-optimized {m(q)) and ^ S{q) /N 
for the ANNNI model at k = 0.7 and qc = 0.1917. The 
diflerence between S{q)/N and {m{q))^ gives xil) (Eg. I20|) . 



4. ORDER AND CRITICAL PARAMETERS 

Structural order parameters help locate phase transi- 
tions, and are particularly important for the study of 
continuous and weakly first-order transitions in models 
studied here. The generalization and application of the 
study of critical and roughening transitions in modulated 
phases is presented in this section. 

4.1. Modulation Order Parameters 

Functions of the Fourier spin density 

N 



i2izqzi 



(17) 



are natural choices for characterizing modulations in lay- 
ered systems. The simplest of them, the generalized mag- 
netization per spin, is defined analogously to the absolute 
magnetization in the Ising model^ 



= {Sq){S-q) 



N 
1 

iV 



\ 



Si cos{qz. 




Sj s'm{qzi) 



A direct use of {m{q)), however, causes problems in long 
simulations, because in principle it averages to zero as the 
lattice drifts (Fig. [3]) . Maximizing the real component of 
Sq with respect to a phase shift in the z direction for each 
configuration before taking the thermal average resolves 
this issue. In practice, we use a straightforward parabolic 
interpolation scheme^. Using the optimized version of 
{m{q)), even in simulations that are too short for the 
system's periodicity to completely diffuse, significantly 
improve the data quality (Fig. [3]). Only quantities based 



on optimized Sq are therefore used in the rest of this work. 
The generalized magnetization decays {m{q)) ~ {T^—T)^ 
with critical exponent /3, but the decay properties are 
not ideal for numerically detecting critical temperatures 
in finite systems. 

The next higher magnetization moment, the z-axial 
static structure factor, is similar to the equivalent liquid- 
state quantity 



NS{q) = {sq~S-q)=N^{m\q)) 



'^SiCOs{qzi) + Si sin{qzi) 



(19) 



where {m?'{q)) is the second moment of the magne- 
tization^iSi'. Both S{q) and its normalized version 
yS{q)/N = ^ {i-n?{q))^^ grow upon cooling and are 
maximal at the wave number qc of the first modulated 
phase below Tc- The monotonically increasing S{q) is, 
however, ill-suited for detecting the PM transition in sim- 
ulations, because, like {m{q)), it does not give a clear 
visual signature of Tc- The generalized susceptibility 

1 



Txiq) = j;^i{sqs-q) - (Sg)(s_g)) 

2 



N{m\q)) -N{m{q)y 



(20) 



i.e., the second cumulant of the magnetization, does not 
suffer from this caveat. It indeed diverges on both sides of 
the transition xilc) ^ \T—Tc\~'' with critical exponent 7, 
as would x(0) in the Ising model, and was used in our pre- 
vious study -° (Fig. [5]). Directly correcting for finite-size 
effects, however, results in a high-sensitivity of the tran- 
sition location to simulation noise. The Binder cumulant 
route is more convenient for detecting Tc, because its 
value at the critical point is straightforwardly insen- 
sitive to scaling the system siz o^^'^°° . For layered phases, 
a generalization of the expression 



Ui{q) = 1 - 



■i{m?{q)Y 



(21) 



in terms of the second and the fourth {m'^{q)) = 
(Sqs'^q) /N'^ q-modulated magnetization moments pro- 
vides the necessary information. 

Because of the anisotropy of the modulated phases, it is 
useful to review how breaking isotropy may affect critical 
properties. In a system of dimensions parallel L|| = Lz 
and perpendicular Lj^ = = Ly to the modulation 
propagation, the correlation length ^ may diverge with 
different critical exponents 



^11 



IT - Tc 



6 



\T-Tc 



(22) 



The critical Binder cumulant 
invariant for a fixed ratio Ln/L 



UiiqcTc) is then 
(Fig. At 



a uniaxial Lifshitz point, such as in the ANNNI model^, 
i^ll ~ 1 ^ ^56,57 ^ ppj. -[-jj^g p]y[ transition, at k > k^, we also 
consider the possibility of anisotropic critical behavior. 
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FIG. 4: Binder cumulant of the ANNNI model at k = 0.8 and 
Qc = 0.2. The curves, which intersect at the critical temper- 
ature Tc — 4.141, monotonically decrease with T. The lim- 
ited validity regime of histogram reweighing results is here re- 
sponsible for the non-monotonicity. (Inset) Finite-size scaling 
analysis of the peak of the derivative of U4{qc) and ln(m(gc)^)- 
The logarithm scales as l/v, which here gives u = 0.66(2). 



Although a direct determination of i'w/va^ is numerically 
difficult, our indirect finite-size study of systems with a 
fixed ratio L\\/L± = 2 shows that 1/^ does not vary at the 
PM transition (Fig. H]) . This observation suggests that 
i'\\/v± ~ 1, i.e., and diverge with the same critical 
exponent = i>± = u at the PM transition. Critical 
anisotropy is thus neglected in the rest of this study. 

Binder cumulants also allow to independently deter- 
mine the critical exponent v using the peak value of the 
derivative of C/4 

In ( -^^-r-] = — lnL + constant. (23) 

A similar relation for the structure factor gives^i^l 
' d\n{m'^{q)) 



In 



dl/T 



1 



■\nL + constant. 



(24) 



The system size L in the scaling relation can be either L|| 
or L± as long as the ratio Ly/Lj, is fixed. For an isotropic 
critical point, as long as the dimensions are rescaled by 
the same factor, the form of the collapse and the critical 
exponents remain unchanged (Fig. |4]). 

Once p and are obtained, the critical exponents (3 
and 7 can more easily be determined through finite-size 
scaling^. The quantities L^^" {m{qc)} and L'^^^^xiQc) 
overlap for different system sizes, when drawn as a func- 
tion of the scaled temperature L^/^iT - Tc)/Tc. The 
heat capacity can also be similarly rescaled, but only 
if C diverges at Tc, as at transitions with Ising univer- 
sality. For a transition with XY-universality, for which 
a = —0.01, C peaks at a finite value in the in- 
finite system size limit. The proper scaling relation is 
then L-^^^iC - C^yo3,io4_ Rgscaling the heat capacity 
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FIG. 5: Finite-size scaling of the magnetization (top), suscep- 
tibility (middle) and heat capacity (bottom) of the ANNNI 
model at K = 0.8 and qc = 0.2, using f = 0.66, 7 — 1.32, 
13 = 0.34, a = -0.01, and = 20. 



curves for such a small a is, however, subject to sizable 
numerical errors^*^'^. The hyperscaling relation 2 — a = 3j/ 
is used instead to determine '"■^ 



4.2. Interfacial Roughening 

In the Ising ferromagnetic regime, even though the cor- 
relation length isotropically diverges at the critical point, 
the interface between two regions of opposite magnetiza- 
tion presents a roughening transition Tr at roughly half 
the critical temperatureiSS^-iSS. Below T/j the interface 
is localized, while above Tr its width diverges logarith- 
mically with surface area. In simulation, an interface is 
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layers whose z coordinates belong to a set /, are grouped 
together in the variance calculation 



FIG. 6: The roughening transition in the ferromagnetic 
regime at k = 0.2 using systems of size x 26 is found 
at Ta = 2.35(5). (insets) Average magnetization profile and 
gradient at T = 3.5. 



created within the bulk by using an antiperiodic bound- 
ary condition along the z directioniS^, and the transition 
can be localized by finite-size analysis (Fig. 15]). Because 
modulated phases intrinsically present a series of inter- 
faces between regions of opposite magnetization, it is also 
interesting to consider whether these interfaces roughen 
or not with temperature. It has been suggested that they 
too should logarithmically diverg o^^i^^" . If that were the 
case, it is possible that interlayer fluctuations at tem- 
peratures between Tr and Tc could participate in phase 
branching and the formation of equilibrium incommen- 
surate structures in the large system limit^ii. A gener- 
alization of the simulation approach is here necessary. 
The variance of the interface position z 



W^^{{z-{z)f) 



(25) 



measures the fluctuations of the interface location, and is 
expected to diverge logarithmically with system size for 
fixed T>Tr 



InL, 



For T < Tr, W'^ should have an even weaker system size 
dependence. In practice, the average in Eq. [25] is taken 
over the normalized magnetization gradient 



{dsxy{z)/dz) 
J (dsxyiz) / dz)dz 

^xy{^ ^" 1) ^xy{^) 
EJSxy(z + l)-S:ry(2)] 



(27) 

(28) 



which serves as weight functionii^. The equilibrium 
profile Sxy{z) is obtained by aligning instantaneous pro- 
files to correct for lattice drift before averaging (Fig. [6]) . 
For the modulated regime, where multiple interfaces are 
present, layers within half a period of the interface i, i.e.. 



zei \zei ) 



(29) 



and the results for the various interfaces are averaged at 
the end. 



5. RESULTS AND DISCUSSION 

Assembling the results from the various observables 
obtained from the computational techniques provides a 
clearer understanding of the equilibrium phase behavior 
of the ANNNI and IC models. In this section, we con- 
centrate on the properties of the modulated regime. 



5.1. Phase Transitions between Layered Phases 

The size of the energy gap between neighboring phases 
with g's commensurate with the simulation box refiects 
the limited and constrained choice of modulations real- 
izable on a finite periodic lattice (Fig. [71). In an infi- 
nite periodic system, where all rational modulations are 
valid but irrational g's are excluded, this gap would be 
infinitely small because rational numbers are dense on 
the real axi o^^'^'^i^^'^ . The smooth and extended energy 
curves for the different modulations are also character- 
istic of strongly metastable phases (Figs. [7| and [8]). The 
high free-energy barriers between layers of differing pe- 
riodicity result in phases that are sufficiently long-lived 
to persist throughout the entire simulation, if the L^Ly 
cross-section is large enough. Different phases can be ob- 
served at a given temperature and frustration, depending 
on how the system is initialized. For smaller cross sec- 
tions, however, the reduced number of spins involved in 
changing the periodicity lowers these transition barriers. 
For a fixed system size, although a longer Lz allows the 
study of more modulated phases, the need to keep these 



(26) 

phases stable limits the maximal aspect ratio Lu : L 



II 

of the simulation lattice. In practice, the selected ratio 
must balance these competing demands. A microscopic 
understanding of the transition mechanism between lay- 
ered phases of different periodicity is still incomplete^, 
but we empirically find that a size ratio of 2:1 is generally 
sufficient. 

The crossing of free-energy curves of neighboring mod- 
ulated phases identifies the transition temperature. Us- 
ing this approach side steps the hysteresis that otherwise 
afflicts annealing approaches, and results in a more ac- 
curate depiction of the modulated regime than had pre- 
viously been obtained^s^i^i^^. For the IC model, for 
instance, the free energy calculations locate the phase 
transitions at temperatures at least 10% lower than re- 
ported in Ref. |43| . where the the system was prepared 
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FIG. 7: (a) Energy and free energy results for the ANNNI model at k = 0.7 for modulations ranging from phase (2) to phase 
(2^3^'*) at melting. The PM transition Tc — 3.988 (vertical dashed line) is obtained from U4. (b) Under the same conditions 
as (a), equilibrium devil's staircase compared with the rescaled soliton result (dashed line), and {m{q)) compared with the 
power- law decay form with /3 = 0.34 obtained from finite-size scaling (solid line). Note that in the low-temperature limit, the 
square profile of phase (2) gives (m{q)) — )■ 2~^^^. (c) Energy and free energy results for the IC model at Q = 0.8 for phases (1), 
(21*) and (21^), and the paramagnetic phase. The PM transition Tc = 1.15(1) is obtained from C. The short vertical dashed 
line indicates the transition temperature between phases (1) and (21^), obtained from simple annealing (see text)^. (d) Under 
the same conditions as (c), devil's staircase and normalized structure factor. Note that in the low-temperature limit, the profile 
of phase (1) gives a/ {m{q)^) — )■ 1. 



in the T = ground state and studied by simulated an- 
nealing. At Q = 0.8, we can even identify a commen- 
surate modulated phase (21^) that was entirely missed 
by the annealing study. The other possibly missed com- 
mensurate phase (21^'^) is, however, unstable here as well, 
presumably because of finite size effects (Fig. [71). Qual- 
itatively similar results are obtained for Q = 0.144 and 
Q = 0.17 (not shown). 

By integrating over frustration at low temperature we 
can also identify the boundary between phases of integer 
periodicity in the IC model. The results at finite temper- 
atures agree very closely with the T = energy derived 
transitions. The location of the free energy cross over be- 
tween phases (1) and (2) (not shown) as well as between 
phases (2) and (3) is only mildly affected by tempera- 
ture (Fig. [8]). The thermal fluctuations produce a similar 
free energy shift of both phases, which leaves the Q loca- 
tion of the transition unchanged. We thus expect similar 
results at other phase {n)-{n + 1) transitions. 



The PM transition is not accurately obtained by di- 
rect free energy comparisons for either systems. For the 
ANNNI model, the continuous transition is best stud- 
ied through the specialized tools of critical phenomena 
(see below). But even for the IC model, the fluctuation- 
induced first-order transition does not lead to sufficiently 
high free energy barriers for noticeably supercooling the 
paramagnetic phase in such a small system. The system 
instead rapidly freezes into a modulated phase below the 
transition and shows only a minimum of hysteresis. As 
a result, the transition identified from the heat capac- 
ity peak by annealing in Ref. is equivalent to what is 
obtained here. A more careful system size dependence 
study would be necessary to refine the transition esti- 
mate. 
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FIG. 8: Energy and free energy curves of the IC model for 
phases (3) and (2) at T = 1.45. The phase boundary from 
free energy calculation (solid lines) agrees with the T = 
mean field prediction (dotted line) and is different from the 
energy inversion (dashed lines)^. 



5.2. Devil's Staircase and Order Parameter 
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The equilibrium wave number obtained from the free 
energy results displays the characteristic devil's stair- 
case2i. The stability regime of a given modulated phase 
stretches over an ever smaller T range upon cooling. For 
the ANNNI model, the predicted truncation of the se- 
quence before reaching the antiphase makes the stair- 
case "harmless" 11^, but the simulated system size is here 
insufficient to distinguish this scenario from the infinite 
"devil's last step" sequence in which no commensurate 
phase is missed^^iii^. The overall shape of the decay 
can, however, be compared with the soliton theory pre- 
diction2i. Though the soliton does not correctly cap- 
ture the PM transition temperature, once T is linearly 
rescaled to make coincide, the agreement is fairly good 
(Fig. El. 

The equilibrium generalized magnetization behaves 
similarly to its g = version in the Ising model. For the 
ANNNI model around Tc, the quantity grows monotoni- 
cally upon cooling. It continuously increases at first, but 
upon reaching the antiphase it jumps discontinuously. In 
the antiphase region the magnetization profile tends to- 
ward a periodic square, whose profile structure is only 
partially captured by a simple sinusoidal function. For 
the IC model, the renormalized structure factor, which 
is indistinguishable from {m{q)) at low temperatures, is 
also not an ideal order parameter. When the system 
changes from phase (21^) to (21''), for instance, the peak 
height actually goes down. Here again, the modulation 
profile is not well captured by a simple sinusoidal func- 
tion. The inclusion of higher order harmonics might bet- 
ter detect growing order upon cooling. 



FIG. 9: Simulation wave number periodicity at melting qc 
for different frustration strengths of the ANNNI (top) and 
the IC (bottom) models superimposed with theoretical pre- 
dictions. The solid line for the IC model captures the ac- 
cessible wave number for Lz — 24. (top inset) The critical 
exponent 7 obtained by finite-size scaling is compared with 
the high-temperature series expansion^, and the field-theory 
predictions for the exponen t^^'^^ . The Lifshitz information is 
taken from Ref. [59. 



5.3. Modulation at Melting 

The periodicity of the modulated phase at the PM 
transition qdn) (or qdQ)) is remarkably insensitive to 
the theoretical approach used for capturing its behavior. 
In the ANNNI model, the agreement between simulation 
results, mean-field theory^, HT series expansion^i*^, 
and the critical scaling near the Lifshitz point 



I A 



(30) 



using either the critical exponent from scries expan- 
sion f3i — 0.5 ± or from renormalization group 
/?/ — 0.514iii, is very good (Fig. [9]). The similarity of 
the RG critical exponent with the mean-field value fur- 
ther suggests that the dependence of microphase period- 
icity on frustration is much easier to capture than the 
transition temperature. The free energy correction due 
to fluctuations is likely similar for neighboring layered 
phases. 

For the IC model, the mean-field prediction for the 
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TABLE I: Critical parameters of the ANNNI model for n < k.l = 0.270(4) obtained by finite-size scaling of cubic systems 
with L — 16, 32, 40, 64, 80 for k = 0.1 and 0.2 and from previous simulations^'^. Ising values are given for reference. The 
uncertainty on Tc and Qc from the HT series expansion results from the Pade approximant method^. At is reported. 

The starred * a results are obtained from the hyperscaling relation = 2 — a (or a + 2/3 + 7 = 2 for k — 0.24). 



K 


Ising 


0.1 


0.15'=^ 


0.2 


24b!> 


0.265'^'' 


0.270''^ 


^ c 


4.512 


4.265(1) 


4.15(2) 


3.987(1) 


3.86(2) 


3.77(2) 


3.7475(5) 


c 


4.51(2) 


4.26(2) 


4.13(2) 


3.98(2) 


3.85(2) 


3.76(2) 


3.75(2) 


V 


0.63 


0.62(1) 


0.61(3) 


0.62(2) 




0.51(4) 


0.33(3)^ 


a 


0.11 


0.14(3)* 0.17(9)* 0.14(6)* 0.28(12)* 0.47(12)* 


0.18(2) 


P 


0.34 


0.31(2) 


0.30(3) 


0.31(3) 


0.23(3) 


0.19(2) 


0.238(5) 


7 


1.24 


1.25(2) 


1.20(6) 


1.23(3) 


1.26(6) 


1.40(6) 


1.36 (3) 



TABLE II: See Table |I] for details. Critical parameters of the ANNNI model for k > k_l obtained by finite-size scaling of 
systems with = 60, 120, 150, 180, and 240 at k = 0.522, = 120, 240, and 360 at k = 0.7, = 60, 100, 120, 200, and 
240 at K = 0.8, = 60, 120, 180, and 240 at k = 2.0, and from previous simulations^!^. XY values are given for reference. 
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(2^3'*) 


(23) 


(2''3) 






0.17(3) 


0.167(4) 


0.18(3) 


0.192(4) 


0.200(4) 


0.233(4) 






0.162(1) 


0.167(1) 


0.180(1) 


0.192(1) 


0.200(1) 


0.232(1) 






0.1667 


0.1705 


0.1816 


0.1919 


0.1994 


0.2301 


rpMC 
^ C 


2.202 


3.6(1) 


3.723(1) 


3.82(3) 


3.988(1) 


4.141(1) 


5.796(1) 


rpHT 

^ c 


2.203ii^ 


3.67(2) 


3.72(2) 


3.81(2) 


3.99(3) 


4.14(1) 


5.79(1) 


V 


0.67 




0.66(2) 




0.67(4) 


0.66(2) 


0.67(3) 


a 


-0.01 




0.02(6)* 




-0.01(9)* 0.02(6)* 


-0.01(9)* 




0.35 




0.35(2) 




0.35(3) 


0.34(3) 


0.36(2) 


7 


1.32 




1.30(4) 




1.32(8) 


1.33(4) 


1.32(6) 



Pade order results in a range of estimates (Tables |T] 
and nil . The values of obtained from finite-size scal- 
ing quantitatively agree with these estimates, and are an 
order of magnitude more precise than both the HT series 
results and previous simulation estimates^-i^ 

The critical exponents from the HT series expansion, 
however, only qualitatively agree with the simulation re- 
sults. Finite HT series can only smoothly approximate 
changes in critical behavior, but the critical exponents 
change discontinuously on both side of the Lifshitz point. 
The HT series results are a continuous approximation of 
that singularity. Field theory arguments suggest, how- 
ever, that the critical exponents should have Ising uni- 
versality below the Lifshitz point, XY above the Lifshitz 
point^-, and uniaxial Lifshitz universality at the Lifshitz 
point^i. The Ising^^ and Lifshitz point^ predictions have 
been previously confirmed by Monte Carlo simulations, 
but above the Lifshitz point the model's behavior is not 
so clear. In particular, the HT series results for 7 at 
large k undershoot the XY exponent value. In the words 
of Ref. ^63, at high k "a puzzling and unexplained feature 
[of the HT series expansion results] is the apparent de- 
crease of 7 to something like the Ising value." Later sim- 
ilar studies did not quite resolve this question, and even 
suggested that a different type of universality might be 
observed beyond ft w 2— . Our earlier simulation results 
did not provide a clear resolution of this issue either, be- 
cause of limited system sizes and insufficient averaging 
in the critical region^^. Simulation of larger systems us- 
ing the multiple histogram method, however, lifts any re- 



continuously changing qc is also within the simulation ac- 
curacy in the layered regime (Fig. ^ , but the relatively 
small lattice size limits quantitatively assessing the the- 
oretical predictions. In the high Q regime, where a Neel- 
paramagnetic-modulated phase triple point is expected, 
our coarse simulation estimate Qn ~ 15.8 clearly dif- 
fers from the mean- field prediction 11.52S. A similarly 
large deviation between the theoretical prediction and 
the direct calculation is also observed at T = 0, where 
Q^/ = 9.549 and 15.33, respectively22i. Both those differ- 
ences can mostly, and possibly completely, be explained 
by the low accuracy of the lattice Fourier transform in the 
large Q limit, where modulated phases of small domains 
fornt2^. Note, however, that the critical nature of the Qn 
point, which depends on the properties of the modulated- 
Neel transition, could also impact its location. If it is a 
bicritical point, fluctuations could result in larger devia- 
tions from the mean- field predictions. A generalization of 
the free energy simulation approach to other modulated 
geometries should be able to resolve this question, but is 
beyond the scope of this work. 



5.4. ANNNI Critical Behavior 

The critical properties of the ANNNI model have been 
extensively studied using high-temperature (HT) series 
expansio n^ 1 ,62 , 70 _ p^j. instance, critical temperatures can 
be estimated by resumming truncated series with Fade 
approximants§2iii^. The arbitrariness of selecting the 
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maining ambiguity. The critical exponent results support 
a XY universality of the transition for all k > kl studied. 
The values of v and 7 agree with each other and with the 
XY values, and are often significantly different from the 
Ising exponents, despite the relatively large error bars 
(Table im and Fig. [9]). The finite-size scaling of C using a 
derived from hyperscaling relations further supports the 
agreement (Fig. [S]). These observations, however, shed 
some doubt on the validity of the predicted transition at 
K « 2. 

The XY universality of the PM transition can be un- 
derstood from the similarity between its two-component 
order parameter and that of the XY mode l^^^i^^° . A 
mean-field picture for the order parameter of the ANNNI 
model (Eq. [T9)) suggests that a spin i in the modulated 
phase can be thought of evolving within a magnetization 
profile of periodicity q formed by all the other spins. A 
change of the average local magnetization at position i is 
equivalent to shifting the phase angle qzi with respect to 
that profile. Note that the isotropic nature of the critical 
PM transition suggests that pairs of spins parallel and 
perpendicular to the z axis are equivalently correlated, 
i.e., the z-axis magnetization profile itself is correlated in 
the X and y directions. In the language of XY model, the 
phase angle is also correlated under translations in the x 
or y directions. 

Why then, one may wonder, do the series expansion re- 
sults not converge to the right 7 value at high k? Exam- 
ining the limit k — > cxd suggests an answer. In that limit 
the next-nearest neighbor interaction dominates and the 
spins decouple into series of intercalated ID Ising anti- 
ferromagnetic chains. That singular limit has ID Ising 
universality for which 7 = 1. The finiteness of the HT 
series expansion thus probably results in a slow decay of 
7 toward unity, as the large k terms in the series dom- 
inate the expansion. In this respect, the series is both 
a high temperature and low k expansion, which further 
restricts its range of validity. 



5.5. ANNNI Roughening Transition 

We first consider the roughening transition of the 
ANNNI model in the ferromagnetic regime (Fig. [6]). 
Though the Tn values extracted from simulations are 
quantitatively different from the series expansion re- 
sultsiSS, similar trends are observed (Figs. [10]). In partic- 
ular, the transition temperature Tr is relatively invariant 
to increases in frustration. The formation of an inter- 
face is not further stabilized by frustration, but rather 
decreases with increasing k. And contrary to the sce- 
nario predicted for other microphase-forming systems, 
the roughening transition does not pass through or near 
the Lifshitz point^^. Instead, the roughening transition 
line on the T-k phase diagram is expected to reach the 
EM phase boundary near k ~ 0.43. Interestingly, a finite- 
temperature intercept suggests that the EM transition 
may be notably different above and below Tr. 



It has also been suggested that a roughening tran- 
sition might be observed for the modulated phases as 
.^g ^iio,i2i _ For the ANNNI model, however, we find no 
indication of interfacial roughening, at least for two sim- 
ple modulated phases: phase (2) at k = 0.8 and phase 
(3) at K = 0.52 (Eig.fTTj). For the latter, in spite of reach- 
ing Tc, the interface location remains clearly defined with 
increasing system size. For any given temperature up to 
the PM transition, the interfacial width of the layers re- 
mains constant upon increasing the system size, but it 
is possible that a divergence can only be observed for 
much larger interfacial areas than what we consider. Yet 
the lattice is here at least an order of magnitude larger 
than the size necessary for detecting roughening in the 
ferromagnetic phase (Fig. |6]). We venture to speculate 
that, at least on a lattice, the persistence length of the 
lamellae might thus be very large and possibly infinite. If 
that were the case, the roughening of the modulated lay- 
ers would then coincide with the PM transition. Further 
simulation and theoretical work are necessary to clarify 
the situation. 



5.6. Phase Diagrams 

Detailed low T series expansion studies of the phase be- 
havior around k = 1/2 conducted by Fisher et. al J^^^^^ 
suggest that a series of "simple phases" of the form (2^ 3) 
spring out from the the multiphase point at T = and 
"mixed phases" generated by combinations of neighbor- 
ing simple phases branch out at T > 0. The tempera- 
tures accessible in simulations are relatively far from the 
regime of validity of this theory and thus, from this point 
of view, it is misleading to compare them directly. It is 
nonetheless interesting to note that the two approaches 
appear to converge for T < 2. 

Various approximate theoretical treatments have been 
used to analyze the ANNNI phase diagram more globally. 
In addition to the traditional mean-field approach''^, an 
effective field^^ and a tensor product variational approach 
(TPVA)^ have more recently been used. These last two 
approaches reasonably capture the external boundaries 
of the modulated regime. The two treatments, how- 
ever, qualitatively disagree on the internal structure of 
that same regime. On the one hand, the effective-field 
method^-, like the mean-field treatment and the soliton 
approximatioEk2i, fills the modulated interior by excep- 
tionally stable bulging simple phases, such as phase (3) 
phase and phase (23) (Fig. [T0|) . On the other hand, 
TPVA predicts rather narrow stability wedges for the 
commensurate phases^^. The simulation results tend to 
favor the second scenario. Though the devil's staircase 
indicates that the rate of wave number change slows on 
approaching Tc, only the antiphase has a broad presence 
in the modulated regime. The stability range of the dif- 
ferent modulations is fairly small, and all of the phases 
commensurate with the periodic box are stable in turn. 
The branching and mixing of the stable low-temperature 
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FIG. 10: Phase boundaries for the ANNNI model obtained from UX (□) and F (0). The Lifshitz point location (•) is taken from 
Ref. [5^. High— and low-temperature^ series expansions as well as TPVA^ results are represented. The stability wedges 
of phases (3) and (23) obtained from simulation are seen to be qualitatively different from the mean-field theory predictions 
(right inset )^. The roughening transition results in the ferromatnetic regime are similar to the series expansion resultsiSS. 




32 X 
64 
128 

256 ' 
320 



5 10 15 20 25 30 35 



0.2 - 




5 10 15 20 25 30 35 



2.5 



3.5 



FIG. 11; Interface width of phase (3) of the ANNNI model 
at ft = 0.52 as a function of T for various system sizes x 
36. No roughening transition is detected within the stability 
regime. The magnetization profile s^^y^z) for the system of 
L = 128 at r = 3.0 is also shown (left inset) along with its 
normalized gradient g(£) (right inset). 



(2-' 3) phases is already complete at the temperatures 
studied here^iii^. In particular, no special stability is 
observed for phases (3) and (23) (Fig.[TU]). For phase (3) 
some bulging is seen, because of the slower rate of change 
of the periodicity near k — 1/2. Simulating larger lat- 
tices, which allow for a more refined q selection, however 
shrinks that phase's footprint. For phase (23), the range 
of stability does increase slightly with k, but the effect is 
probably due to the finiteness of the lattice. In any case. 



the increase is much less pronounced than the bulging 
scenarios predicl—iS (Fig. ITU)). 

For the IC model, qualitatively similar branching is 
expected in the devil's flower region, which is found be- 
tween phases of integer periodicity. The systems simu- 
lated here are, however, too small to examine this issue 
critically. Except for the caveats presented in the previ- 
ous sections, our results mostly agree with the simulation 
results of Ref. IS 



6. CONCLUSION 

Our simulation study has clarified the structure and 
transition properties of the modulated regime of the 
ANNNI and the IC models. Previous theoretical treat- 
ments had sometimes been insufficient, particularly con- 
cerning the stability regime of the various modulated 
phases, the critical nature of the ANNNI PM transition, 
and the role of roughening. In the last case, no clear 
conclusion can be drawn, but the results suggest that 
the phenomenon is at least a lot less pronounced than in 
the Ising model, which may give hope of experimentally 
forming microphase patterns on much larger scales than 
previously thought. From a theoretical perspective, it 
is also interesting to highlight, however, that mean-field 
theory is particularly adept at predicting the periodicity 
of modulated phases at the PM transition. This observa- 
tion may explain why the approach has been so successful 
at describing order in other microphase-forming systems, 
such as diblock copolymerai^ii2^. 

In addition to lamellar phases, modulated assemblies 
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can exhibit a variety of other symmetries. They can 
also be observed off lattice. Generalizing the approach 
to continuous space and to other order types would thus 
greatly benefit the study of more complex microphase- 
forming systems. For the IC model, it could for in- 
stance help determine the nature of the modulated-Neel 
transition and other properties of the high Q regime, 
which we only briefly explored. Completing the simu- 
lation tool set would also pave the way for studies of the 
non-equilibrium microphase assembly, where most of the 



materials challenges lie. 
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